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For the first time weak turbulent theory was demonstrated for the surface gravity waves. Direct numerical 
simulation of the dynamical equations shows Kolmogorov turbulent spectra as predicted by analytical analysis 
1^ from kinetic equation. 
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^ ' In this Letter we study numerically the steady 
^ iKolmogorov spectra for spatially homogeneous gravity 

^ I I waves. According to the theory of weak turbulence 
the main physical process here is the stationary en- 

_ ^ ergy flow to the small scales, where the energy dissi- 
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pates 13 . This flow is described by kinetic equation 
'which has power-like solutions - Kolmogorov spectra. 
This straightforward picture takes place experimentally 
land numerically for different physical situations. For 
■capillary waves it was observed on the surface of liq- 
uid hydrogen 0, The numerical simulation of this 
process was performed in [HI. In nonlinear fiber optics 
.these spectra were demonstrated in numerical simula- 
tion ini-There are many other results jZl El El HIB CI] • 
One of the most interesting applications of the weak 
turbulence theory is the surface gravity waves. From 
the pioneering article by Toba ^21 to the most recent 
^ .observations many experimentalists get the spectra 
. ^ 'predicted by the weak turbulence theory. But these ex- 
^_^,periments cannot be treated as a complete confirmation 
(~| jbecause the Zakharov-Filonenko spectrum is isotropic, 
Q^iwhile observed spectra are essentially anisotropic. It is 
^ worth to say that the wave kinetic equation, which is 
! the keystone of this theory, was derived under several 
I^S assumptions. Namely, it was assumed, that the phases 
of all interacting waves are random and are in state of 
chaotic motion. The validity of this proposition is not 
clear a priori. The direct numerical simulation of non- 
linear dynamical equations can give us a confirmation 
is this assumption valid or not. But for particular case 
of gravity surface waves the numerical confirmation was 
absent in spite of significant efforts were applied. The 
only successful attempt in this direction was the sim- 
ulation of freely decaying waves jEj. The reason for 
that for our opinion was concerned with a choice of 
numerical scheme parameters. Namely, the numerical 



simulation is very sensitive to the width of resonance 
of four-waves interaction. It must be wide enough to 
provide resonance on the discrete grid, as it was studied 
in 15 for decay of the monochromatic capillary wave. 
From the other hand it has to be not too wide (due to 
nonlinear frequency shift) when the weak turbulent con- 
ditions fail. We have spent significant efforts to secure 
the right choice of numerical parameters. As a result we 
have obtained the first evidence of the weak turbulent 
Kolmogorov spectrum for energy flow for surface grav- 
ity waves. The numerical simulation was surprisingly 
time consuming (in comparison to capillary waves tur- 
bulence) , but finally we clearly get spectrum for surface 
elevation ^ 



which is in the agreement with real experiments |12II13| . 

Theoretical background. — Let us consider the po- 
tential flow of an ideal incompressible fluid of infinite 
depth and with a free surface. We use standard nota- 
tions for velocity potential 0(r, z,t),r = {x,y);v = V(/> 
and surface elevation ri{r,t). Fluid flow is irrotational 
A0 — 0. The total energy of the system can be repre- 
sented in the following form 

H = T + U, 



T=- I d\ I {^4>fdz, 



(2) 



(3) 
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where g - is the gravity acceleration. It was shown 
that under these assumptions the fluid is a Hamiltonian 
system 

d7j _^SH di}j _ SH 
dt Sip ' dt 6r] ' 



(4) 
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where ip — (l){r,ii{r,t),t) is a velocity potential on the 
surface of the fluid. In order to calculate the value of 
^ we have to solve the Laplas equation in the domain 
with varying surface 77. This problem is difficult. One 
can simplify the situation, using the expansion of the 
Hamiltonian in powers of " steepness" 

in-'' 
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For gravity waves it is enough to take into account terms 
up to the fourth order. Here k is the linear operator 
corresponding to multiplying of Fourier harmonics by 
modulus of the wavenumber k. In this case dynamical 
equations ^ acquire the following form 



j'y = k^ — (V(?7V'0)) — k[rjk^ 
+k{r]k['qk^]) + \l\[rf'ki} 



ip = -gr] 



(6) 



[kxp]k[T]kij] - [rikTp]A'ilj + + 



Here is some artificial damping term used to provide 
dissipation at small scales; Fr is a pumping term corre- 
sponding to external force (having in mind wind blow, 
for example). Let us introduce Fourier transform 

1 r . 1 
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With these variables the Hamiltonian lO acquires the 
following form 



H = Ha 



Ho = 


A-K 


Hi = 
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H2 = 
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-?-kik2'0kiV'k2?7k3X 

k3)dki(ik2dk3, 

A^kikskakiV'ki V'k2'7k3'7k4 X 

ka + k4)dkidk2dk3C?k4, 



(7) 



Here 

ikika = (kik2) + \ki 
-M"kik2k3k4 = |ki||k2| 

+ |k2+k3| + |k2 
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ki +k3| + |ki +k4| 
-|ki|-|k2|]. 



(8) 



It is convenient to introduce the canonical variables 
as shown below 



l^k k 

2fc^'' + V^^''' 



(9) 



where 

ujk = \/gk, (10) 

this is the dispersion relation for the case of infinite 
depth. The similar formulas can be derived in the case 
of finite depth ^|. With these variables the equations 
Q take the following form 

.SH 
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The dispersion relation (|10() is of the " non-decay type" 
and the equations 



LUki = U!k2 + ^k:, , ki 



(12) 



have no real solution. It means that in the limit of 
small nonlinearity, the cubic terms in the Hamiltonian 
can be excluded by a proper canonical transformation 
a(k, t) — > &(k, t) JHl- The formula of this transforma- 
tion is rather bulky and well known |2I ^| , so let us 
omit the details here. 

For statistical description of a stochastic wave field 
one can use a pair correlation function 



< akflk' nk5{h — k'). 



(13) 



The rik is measurable quantity, connected directly with 
observable correlation functions. For instance, from @ 
one can get 

4. -<hkP>=i — (nfe-Hn^fc). (14) 
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In the case of gravity waves it is convenient to use an- 
other correlation function 



<6k6k' >^NkSik-k'). 



(15) 



The function Nk cannot be measured directly. The rela- 
tion connecting and Nk is rather complex in the case 
of fluid of finite depth. But in the case of deep water it 
becomes very simple |17j 

rik - Nk 



rik 



(16) 



where /i = (ka)'^, here a is a characteristic elevation 
of the free surface. In the case of the weak turbulence 
II « 1. The correlation fimction Nk obey the kinetic 
equation P 

dNk 



dt 



Here 



= st(iV,7V,iV) + /p(fc)-/rf(fc), (17) 



(18) 



st{N,N,N)=A^ J |Tk,ki,k2,k3rx 
x{Nk,Nk,Nk, + NkNk.Nk, - NkNk,Nk,~ 

~NkNk,Nk,)6ik + ki - k2 - k3)dkidk2dk3. 
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The complete form of matrix element T'k.ki,k2,k3 can be 
found in many sources QIEIEI- Function fp(k) in H17|) 
corresponds to wave pumping due to wind blow for ex- 
ample. Usually it is located on long scales. Function 
fd{k) represents the absorption of waves due to viscos- 
ity and wave-breaking. None of this functions are known 
to a sufBcient degree. 

Let us consider stationary solutions of the equation 
ifTTjl assuming that 

• The medium is isotropic with respect to rotations; 

• Dispersion relation is a power-like function lu — 
ak"; 

• 7k.ki,k2,k3 is a homogeneous function: 

^ek.eki ,ek2 ,€k3 — ^^^k.ki ,k2 ,k3 . 

Under this assumptions one can get Kolmogorov solu- 
tions [IS] 



(2) . CQ'/^k-'-^ 



(19) 



Here c? is a spatial dimension (c? = 2 in our case). 
The first one is a Kolmogorov spectrum, correspond- 
ing to a constant flux of energy P to the region of small 
scales (direct cascade of energy). The second one is Kol- 
mogorov spectrum, describing inverse cascade of wave 
action to large scales, and Q is a flux of action. In both 
cases Ci and C2 are dimensionless "Kolmogorov's con- 
stants" . 

In the case of deep water uj = \/gk and, apparently, 
/3 = 3. It is known since 2] that on deep water 



CiP'/^k-\ 



In the same way |19| for second spectrum 



(20) 



(21) 



In this Letter we will explore the first spectrum (en- 
ergy cascade). Using ifT^ one can get 



Ik 



(22) 



Numerical Simulation — Dynamical equations © 
are very hard for analytical analysis. One of the main 
obstacles is the fc-operator which is nonlocal. However, 
using Fourier technique practically makes no difference 
between derivative and k. The numerical simulation of 
the system is based upon consequent application of fast 
Fourier transform algorithm. The details of this numer- 
ical scheme will be published separately. 



For numerical integration of ^ we used the func- 
tions F and D deflned in Fourier space 



(23) 



= ^j.^ (fc-fcpl)(fcp2-fe) , 

{kp2-kpiY 

Ik = -li-k < kpi, 

Ik = -l2{k - kdf, k > kd- 



Here Rk{t) is the uniformly distributed random num- 
ber in the interval (0, 27r). We have solved system of 
equations ^ in the periodic domain 27r x 27r (the wave- 
numbers kx and ky are integers in this case). The size 
of the grid was chosen 256 x 256 points. Gravity accel- 
eration g = 1. Parameters of the damping and pumping 
were the following: kpi — 5, kp2 — 10, kd ~ 64. Thus 
the inertial interval is about half of decade. 

During the simulations we paid special attention 
to the problems which could "damage" the calcula- 
tions. First of all, the "bottle neck" phenomenon at the 
boundary between inertial interval and dissipation re- 
gion. This effect is very fast, but can be effectively sup- 
pressed by proper choice of damping value 72 in the case 
of moderate pumping values Fq. The second problem is 
the accumulation of " condensate" in low wave numbers. 
This mechanism for the case of capillary waves was ex- 
amined in details in ^HI- This obstacle can be over- 
come by simple adaptive damping scheme in the small 
wave numbers. After some time system reaches the sta- 
tionary state, where the equilibrium between pumping 
and damping takes place. Important parameter in this 
state is the ratio of nonlinear energy to the linear one 
{H, + H2)/Ho. 

For example, in the case of Fg = 2 x 10~^,7i = 
1 X 10"'^, 72 = 400 the level of nonlinearity was equal 
to {Hi + H2)/Hq ~ 2 X 10~^. The Hamiltonian as a 
function of time is shown in Fig. ^ 

The surface elevation correlator function appears to 
be power-like in the essential part of inertial interval, 
where the influence of pumping and damping was small. 
The correlator is shown in Fig. [2| 

One can try to estimate the exponent of the spec- 
trum. It is worth to say that an alternative spectrum 
was proposed earlier by Phillips [201 • That power- like 
spectrum is due to wave breaking mechanism and gives 
us a surface elevation correlator as Ik ^ Compen- 
sated spectra are shown in the Fig. |3| It seems to be 
an evidence, that the Kolmogorov spectrum predicted 
by weak turbulence theory better fit the results of the 
numerical experiment. 
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Fig.l. Hamiltonian as a function of time. 
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Fig.2. The logarithm of the correlator function of sur- 
face elevation as a function of logarithm of the wave 
number. 



The inertial interval was rather narrow (half a 
decade). But the obtained results allow us to conclude, 
that accuracy of experiment was good enough under the 
time constraints of simulation (we get the steady state 
after 20-30 h using available hardware, and we need sev- 
eral days to average |?7fcp function). The simulation on 
larger grid (512 x 512, for example) can make the ac- 
curacy better. But even these results give us a clear 
qualitative picture. 
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Fig. 3. Compensated correlators in inertial interval for 
different values of the compensation power; z — 3.5 
solid line (weak turbulence theory), z = 4.0 dashed line 
(Phillips theory). 
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